%pylab inline
import matplotlib.pyplot as plt

from scipy.optimize import fsolve

from ipywidgets import *

Populating the interactive namespace from numpy and matplotlib

Relativistic kinematics (Relativisztikus kinematika)¶

Decay of a Particle to Two Daughters: $A \to B+C$¶

Eqs. for decay:¶

\begin{eqnarray} E_A &=& E_B+E_C, \\ p_A &=& p_B+p_C, \\ E_A^2/c^2 -p_A^2 &=& M_A^2 c^2,\\ E_B^2/c^2 -p_B^2 &=& M_B^2 c^2,\\ E_C^2/c^2 -p_C^2 &=& M_C^2 c^2. \end{eqnarray}

Input: $p_A,M_A,M_B,M_C$¶

Output: $E_A,E_B,E_C, p_B,p_C$¶

The solution (in units of $c=1$):¶

\begin{eqnarray} E_A &=& \sqrt{p_A^2+M_A^2}, \\ E_B &=& \frac{1}{2 M_A^2}\left[-Q p_A +\left(M_A^2+M_B^2-M_C^2\right) E_A \right], \\ E_C &=& \frac{1}{2 M_A^2}\left[Q p_A +\left(M_A^2-M_B^2+M_C^2\right) E_A \right], \\ p_B &=& \frac{1}{2 M_A^2}\left[-Q E_A +\left(M_A^2+M_B^2-M_C^2\right) p_A \right],\\ p_C &=& \frac{1}{2 M_A^2}\left[Q E_A +\left(M_A^2-M_B^2+M_C^2\right) p_A \right], \textrm{where} \\ Q &=& \sqrt{\left[M_A^2-{\left(M_B+M_C\right)}^2\right]\left[M_A^2-{\left(M_B-M_C\right)}^2\right]}, \end{eqnarray}

and the velocities are calculated from $p= E V/c^2$.¶

def decay(be):
    
    '''
    Analitic solution for the above equations
    
    input: pA,MA,MB,MC
    output: EA,EB,EC,pB,pC
    
    '''
    
    pA,MA,MB,MC = be 
     
    EA=sqrt(pA**2+MA**2)
    #Q= sqrt(MA**4 + (MB**2 - MC**2)**2 - 2*MA**2 *(MB**2 + MC**2) )
    Q= sqrt((MA**2 - (MB - MC)**2)*(MA**2 - (MB + MC)**2))
    EB=1/2/MA**2*(-Q*pA+(MA**2+MB**2-MC**2)*EA)
    EC=1/2/MA**2*(Q*pA+(MA**2-MB**2+MC**2)*EA)
    pB=1/2/MA**2*((MA**2+MB**2-MC**2)*pA-Q*EA)
    pC=1/2/MA**2*((MA**2-MB**2+MC**2)*pA+Q*EA)
    
    return(EA,EB,EC,pB,pC)

def func_mass_shell(MA,chimax,szin):
    '''
    plotting the mass shell for p and E
    
    '''
    chi=linspace(-chimax,chimax,100)
    plot(MA*sinh(chi),MA*cosh(chi),'-',color=szin);

def rajz_mass_shell(pA,paramM):  
    '''
    plotting the vector (p,E) for the elements A, B and C
    '''
    
    MA,MB,MC = paramM
    
    figsize(10,6)
    #ax=subplot(aspect='equal')
    ax=subplot(111,aspect='equal')

    chimax=[0.7,1,2.2]
    func_mass_shell(MA,chimax[0],'b')
    text(1.01*MA*sinh(chimax[0]),1.01*MA*cosh(chimax[0]),'A',color='b',fontsize=20);
    
    func_mass_shell(MB,chimax[1],'r')
    text(1.01*MB*sinh(chimax[1]),1.01*MB*cosh(chimax[1]),'B',color='r',fontsize=20);
    
    func_mass_shell(MC,chimax[2],'g')
    text(1.01*MC*sinh(chimax[2]),1.01*MC*cosh(chimax[2]),'C',color='g',fontsize=20);
    
    fact=1.2
    
    plot([0,fact*MA],[0,0],'-k')
    text(1.05*fact*MA,0,'$p$',fontsize=20)
    
    plot([0,0],[0,fact*MA],'-k')
    text(0,1.05*fact*MA,'$E$',fontsize=20)
    #plot([0,fact*MA],[0,fact*MA],'--k')
    
    param = [pA]+paramM   # adding two arrays  
    EA,EB,EC,pB,pC = decay(param)   # res --> EA,EB,EC,pB,pC
    
    print('p_A = ',pA,'E_A = ',EA,'V_A = ',pA/EA,'M_A = ',MA,'\n')
    print('p_B = ',pB,'E_B = ',EB,'V_B = ',pB/EB,'M_B = ',MB,'\n')
    print('p_C = ',pC,'E_C = ',EC,'V_C = ',pC/EC,'M_C = ',MC,'\n')

    print('mass defect =',param[1]-param[2]-param[3])
    
    plot([0,pA],[0,EA],'-b')
    plot([0,pB],[0,EB],'-r')
    plot([0,pC],[0,EC],'-g')
    
    plot([pC,pB+pC],[EC,EB+EC],'--r')
    plot([pB,pB+pC],[EB,EB+EC],'--g')
    
    title(r'$A \to B + C$',fontsize =20)
    
    grid();

Table for rest mass of mesons ¶

$K^{*-} \to K^{-} + \pi^0$¶

param =  [0.0,0.8917,0.4937,0.1350]  # pA,MA,MB,MC
pA,MA,MB,MC =param

EA,EB,EC,pB,pC = decay(param)   # res --> EA,EB,EC,pB,pC

print('p_A = ',pA,'E_A = ',EA,'V_A = ',pA/EA,'M_A = ',MA,'\n')
print('p_B = ',pB,'E_B = ',EB,'V_B = ',pB/EB,'M_B = ',MB,'\n')
print('p_C = ',pC,'E_C = ',EC,'V_C = ',pC/EC,'M_C = ',MC,'\n')
    
print('mass defect =',param[1]-param[2]-param[3])

p_A =  0.0 E_A =  0.8917 V_A =  0.0 M_A =  0.8917 

p_B =  -0.2894650465975274 E_B =  0.5723021083323988 V_B =  -0.5057906346719296 M_B =  0.4937 

p_C =  0.2894650465975274 E_C =  0.3193978916676012 V_C =  0.9062835233075833 M_C =  0.135 

mass defect = 0.263

Neutron Decay: $n \to p + e^{-} + \bar{\nu}$¶

param =  [0,939.5656, 938.2723, 0.5110]  # pn,Mn,Mp,Me
pA,MA,MB,MC =param

EA,EB,EC,pB,pC = decay(param)   # res --> EA,EB,EC,pB,pC

print('p_A = ',pA,'E_A = ',EA,'V_A = ',pA/EA,'M_A = ',MA,'\n')
print('p_B = ',pB,'E_B = ',EB,'V_B = ',pB/EB,'M_B = ',MB,'\n')
print('p_C = ',pC,'E_C = ',EC,'V_C = ',pC/EC,'M_C = ',MC,'\n')
    
print('mass defect =',param[1]-param[2]-param[3])

p_A =  0 E_A =  939.5656 V_A =  0.0 M_A =  939.5656 

p_B =  -1.187249568211231 E_B =  938.2730511470675 V_B =  -0.001265356142073762 M_B =  938.2723 

p_C =  1.187249568211231 E_C =  1.2925488529327245 V_C =  0.9185336132691813 M_C =  0.511 

mass defect = 0.7823000000000447

$K^{0} \to \pi^+ + \pi^-$¶

param =  [667,500,140,140] # pA,MA,MB,MC
pA,MA,MB,MC =param

EA,EB,EC,pB,pC = decay(param)   # res --> EA,EB,EC,pB,pC

print('p_A = ',pA,'E_A = ',EA,'V_A = ',pA/EA,'M_A = ',MA,'\n')
print('p_B = ',pB,'E_B = ',EB,'V_B = ',pB/EB,'M_B = ',MB,'\n')
print('p_C = ',pC,'E_C = ',EC,'V_C = ',pC/EC,'M_C = ',MC,'\n')
    
print('mass defect =',param[1]-param[2]-param[3])

p_A =  667 E_A =  833.6000239923221 V_A =  0.8001439309053372 M_A =  500 

p_B =  -11.81572857314217 E_B =  140.4977275322066 V_B =  -0.08409907249520193 M_B =  140 

p_C =  678.8157285731421 E_C =  693.1022964601154 V_C =  0.9793875046152073 M_C =  140 

mass defect = 220

pA=0.27
paramM = [0.8917,0.4937,0.1350]   # MA,MB,MC
rajz_mass_shell(pA,paramM)

p_A =  0.27 E_A =  0.9316806802762414 V_A =  0.28979886104319086 M_A =  0.8917 

p_B =  -0.12915489770150335 E_B =  0.5103142929610006 V_B =  -0.2530889286915851 M_B =  0.4937 

p_C =  0.39915489770150336 E_C =  0.421366387315241 V_C =  0.9472869923126537 M_C =  0.135 

mass defect = 0.263

paramM = [0.8917,0.4937,0.1350]   # MA,MB,MC

@interact(pA=(-1,1,0.01))

def play(pA=0.2):
    
    rajz_mass_shell(pA,paramM)

MA,MB,MC = [0.8917,0.4937,0.1350]
Q= sqrt((MA**2 - (MB - MC)**2)*(MA**2 - (MB + MC)**2))
EA=0.9138538668736923 
#pA=EA*Q/(MA**2+MB**2-MC**2)
VA=Q/(MA**2+MB**2-MC**2)
VA

0.5057906346719295